%
% Example generation
%

%
% signal to noise ratios
SNRx=0;
SNRy=0;
% Coupling parameter
beta=.5;
%
%
N=1024; % number of data points
e=rand(2,N);
x=zeros(N+2,1);
y=x;
% Dynamical coupling
for k=3:N+2
    x(k)=1.5806*x(k-1)-.7225*x(k-2)+e(1,k-2);
    y(k)=.5*y(k-1)+e(2,k-2)+ beta*x(k-1);
end
x=x(3:N+2);
y=y(3:N+2);
% square root powers
sx=std(x);
sy=std(y);
%
Ax=sx*10^(-SNRx/20);
Ay=sx*10^(-SNRx/20);
%
xu=Ax*cos(2*pi*60*(0:N-1)/250);
yu=Ay*cos(2*pi*60*(0:N-1)/250);
% final signal
xt=x+xu';
yt=y+yu';
%
xt=xt-mean(xt);
yt=yt-mean(yt);
% notch filter parameters for BW=.5 Hz - 250 Hz sampling for 60 Hz.
num=[ 0.9978   -0.1253    0.9978];
den=[1.0000   -0.1253    0.9956];
%
xf=filter(num,den,xt);
yf=filter(num,den,yt);
%
figure
plot(x)
hold
plot(y,'r')
hold
figure
plot(xt)
hold
plot(yt,'r')
hold
figure
plot(xf)
hold
plot(yf,'r')
hold


